// project home: cp3-cg.googlecode.com

#include "cp31.h"
#include <Rtypes.h>
#include <iostream>
#include <cstdlib>
#include <random>
#include <vector>
#include <time.h> 
#include "cg.h"
#include "myvectormath.h"
using namespace std;


Int_t main ( Int_t argc, const char* argv[] ){
  
  // Prints each argument on the command line.
  for( int i = 0; i < argc; i++ )
  {
    cout << "arg " << i << ": " << argv[i] << endl;
  }
  
  
  // create matrix
  // get dimension of matrix
  Int_t fMatDim;
  int debug;
  if (argc==1){cout << "Call with \"cp31 MatDim debugv\" \n MatDim = Dimension of the Matrix \n debugv = debug verbosity" << endl;}
  if (argc>1){fMatDim=atoi(argv[1]);}else{fMatDim=5;}
  if (argc>2){debug=atoi(argv[2]);}else{debug=0;}
  
  if (debug){cout << "Matrix Dimension: " << fMatDim << endl;}
  
  
  
  // generators for non-zero elements of matrix. for now only diagonal matrices are created and so the position is fixed.
  
  default_random_engine generator(time(NULL));
  uniform_real_distribution<Double_t> distributiond(0.0,1.0); // to generate the values
  // uniform_int_distribution<int> distributioni(1,fMatDim); // to generate the position of the nonzero elements
  
  // The matrix is stored with compressed row storage, compare: http://www.ma.utexas.edu/CNA/ITPACK/manuals/user2c/node2.html
  // Matrix A holds the values of the nonzero elements, Matrix JA holds the column of these elements and Matrix IA stores when the next row starts
  
  /*
  // fMatDim of entries are != 0; creates diagonal matrix with real elements greater or equal zero 
  vector<Double_t> A(fMatDim);
  vector<Double_t> JA(fMatDim);
  vector<Double_t> IA(fMatDim+1);
  for( int i = 0; i < fMatDim; i++) {
    A.at(i)= distributiond(generator);
    JA.at(i)=i;
    IA.at(i)=i;
  }
  IA.at(fMatDim)=fMatDim;
  */
  
  // creating matrix with ca. 2*N elements at random positions
  vector<Double_t> B;
  vector<Int_t> JB;
  vector<Int_t> IB;
  Double_t probelem=(Double_t)2/(Double_t)fMatDim;
  int testcount=1;
  for (int i=0; i<fMatDim; i++){
    IB.push_back(B.size());
    for (int j=0; j<fMatDim; j++){
//     if(distributiond(generator)<probelem){
      if(true){
	B.push_back(distributiond(generator));
// 	B.push_back(testcount);
	JB.push_back(j);
	testcount++;
      }
    }
  }  
  IB.push_back(B.size());
  
  if (debug){cout << "creating the first matrix done" << endl;}
  
  vector<Double_t> rowin(fMatDim,0);
  vector<Double_t> rowout(fMatDim);
  vector<Double_t> A;
  vector<Int_t> JA;
  vector<Int_t> IA;
  for (int row=0;row<fMatDim;row++){
    rowin.assign(fMatDim
    ,0);
    if (debug>=3){cout << "row:" << row << endl;}
    IA.push_back(A.size());
    for (int col=IB.at(row); col<=IB.at(row+1)-1; col++){
      if (debug>=3){cout << "col" << col << endl;}
      rowin.at(JB.at(col))=B.at(col);
    }
    rowout=matrixvectorproduct(&B,&JB,&IB,&rowin);
    for (int colout=0; colout<fMatDim; colout++){
      if (rowout.at(colout)!=0){
	A.push_back(rowout.at(colout));
	JA.push_back(colout);
      }
    }
  }
  
  IA.push_back(A.size());
  

  
  
  
  if (debug){cout << "creating the matrix done" << endl;}
  
  
  //debug the creating of the matrix
  if (debug>=2){
    // print the vectors and number of elements in the matrix
    cout << "testcount: " << testcount-1 << endl;
    cout << "A: ";
    for (auto it = A.cbegin(); it != A.cend(); ++it)
      std::cout << ' ' << *it;
    std::cout << '\n';
    
    cout << "JA: ";
    for (auto it = JA.cbegin(); it != JA.cend(); ++it)
      std::cout << ' ' << *it;
    std::cout << '\n';
    
    cout << "IA: ";
    for (auto it = IA.cbegin(); it != IA.cend(); ++it)
      std::cout << ' ' << *it;
    std::cout << '\n';
    
    
    
    // first: transform the three vectors in a series of 3-tupels with (row, column, value)
    cout << "transforming the matrix into 3-tupels: " << endl;
    vector<vector<Double_t> > Atupel;
    vector<Double_t> elemtupel(3);
    vector<Double_t>::iterator itA = A.begin();
    vector<Int_t>::iterator itJA = JA.begin();
    for (int row=0; row<fMatDim; row++){
      if (debug>=3){cout << "row=" << row << endl;}
      if (debug>=3){cout << "IA.at(row)=" << IA.at(row) << endl;}
      if (debug>=3){cout << "IA.at(row+1)=" << IA.at(row+1) << endl;}
      for (int col=IA.at(row); col<=IA.at(row+1)-1;col++){
	if (debug>=3){cout << "col=" << col << endl;}
	elemtupel.at(0)=row;
	elemtupel.at(1)=*itJA;
	elemtupel.at(2)=*itA;
	itA++;
	itJA++;
	Atupel.push_back(elemtupel);
      }
    }
    
    // print the 3-tuples
    cout << "Atuples: ";
    for (auto it = Atupel.cbegin(); it != Atupel.cend(); ++it)
      std::cout << "(" << it->at(0) << "," << it->at(1) << "," << it->at(2) << ") ";
    std::cout << '\n';
    
    
    
    // second: print out that matrix
    cout << "printing the matrix: " << endl;
    vector<vector<Double_t> >::iterator itAtupel = Atupel.begin();
    for (int i=0; i<fMatDim; i++){
      for (int j=0; j<fMatDim; j++){
	if( itAtupel!= Atupel.end() ){
	  if (i==itAtupel->at(0)&&j==itAtupel->at(1)){
	    cout << itAtupel->at(2) << " ";
	    itAtupel++;
	  }else{
	    cout << 0 << " ";
	  }
	}else{
	  cout << 0 << " ";
	}
      }
      cout << endl;
    }

  }
  
  
  
  // generate the vector b
  vector<Double_t> b(fMatDim);
  for( int i = 0; i < fMatDim; i++) {
     b.at(i)= 2*distributiond(generator)-1;
//    b.at(i)= 2;
  }
  
  if (debug){cout << "creating b done" << endl;}
  
  vector<Double_t> testv(fMatDim);
  
  testv=(matrixvectorproduct(&A,&JA,&IA,&b,debug));
  
  if (debug){cout << "product between matrix and vector done" << endl;}
  
  if (debug){
    cout << "b: ";
    for (int i=0;i<fMatDim;i++){
      cout << b.at(i) << " ";
    }
    cout << endl;
    cout << "testv: ";
    for (int i=0;i<fMatDim;i++){
      cout << testv.at(i) << " ";
    }
    cout << endl;
  }
  
  
  vector<Double_t> x(fMatDim,0);
  cout << cg(&A, &JA, &IA, fMatDim, &x, &b, "fRoutine", pow(10,-8), fMatDim*2, true, debug);
  //cg.h(&A, &JA, &IA,fMatDim, &x, &b
    
  return 0;
}




bool matelemexists(vector<Int_t> *JA, vector<Int_t> *IA, Int_t row, Int_t col){
  for (int i=IA->at(row); i<=IA->at(row+1)-1;i++){
    if (JA->at(i)==col){ return true; }
  }
  return false;
}






